Mixed Models and GxE
A tutorial on how to run asses genotype by environment interactions in R
Introduction
- GxE: https://en.wikipedia.org/wiki/Gene%E2%80%93environment_interaction -Mixed Models https://en.wikipedia.org/wiki/Mixed_model
https://www.azandisresearch.com/2022/12/31/visualize-mixed-effect-regressions-in-r-with-ggplot2/
https://nariyoo.com/stata-introduction-to-multilevel-modeling-xtmixed/
https://bookdown.org/steve_midway/DAR/random-effects.html
# devtools::install_github("derekmichaelwright/agData")
library(agData)
library(lme4) # For mixed-effects models
# Prep data
myCaption <- "derekmichaelwright.github.io/dblogr/ | Data: AGILE"Data
- global-energy-substitution.csv
- Wright, et al. Understanding photothermal interactions can help expand production range and increase genetic diversity of lentil (Lens culinaris Medik.). Plants, People, Planet. (2020) 00: 1-11. doi.org/10.1002/ppp3.10158
GxE
GxE and Plotting Functions
myGxE <- function(xx, myTraits = c("DTE","DTF","DTS","DTM","VEG","REP")) {
# Initialize an empty list to store variance component data
var_list <- list()
blup_list <- list() #store blups values
# Loop over each trait
for(i in myTraits) {
# Define the model formula with ENV as a fixed factor
formula <- as.formula(paste(i, "~ Expt + (1 | Name) + (1 | Name:Expt) + (1 | Expt:Rep)"))
# Fit the mixed-effects model
model <- lmer(formula, data = xx) # Replace with your actual data
# Extract variance components for random effects
var_i <- VarCorr(model) %>% as.data.frame(vc)
# Calculate total phenotypic variance
var_total <- sum(var_i$vcov)
# Calculate proportions as percentages
var_i$Proportion <- (var_i$vcov / var_total) * 100
#
var_i <- var_i %>% mutate(Trait = i) %>%
select(Component=grp, Trait, Proportion, vcov, sdcor)
# Store results in the list
var_list[[i]] <- var_i
}
# Combine all trait data into one data frame
var_table <- do.call(rbind, var_list) %>%
mutate(Trait = factor(Trait, levels = myTraits),
Component = factor(Component,
levels = c("Residual", "Expt:Rep", "Name:Expt", "Name"),
labels = c("Residual + Env", "Rep:Env", "G x E", "Genetic")))
# Output
var_table
}gg_GxE <- function(xx) {
myColors <- c("darkred", "purple4", "steelblue", "darkgreen")
# Create the stacked bar plot
mp <- ggplot(xx, aes(x = Trait, y = Proportion, fill = Component)) +
geom_col(position = "stack", color = "black", alpha = 0.7) +
scale_fill_manual(values = myColors) +
theme_agData_col() +
labs(title = "Proportion of Phenotypic Variance by Trait",
y = "Proportion of Total Variance (%)",
caption = myCaption)
}All Environments
# Save variance table
var_table <- myGxE(dd)
write.csv(var_table, "var_All.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("variance_All.png", mp, width = 6, height = 4)Temperate
# Save variance table
var_table <- myGxE(dd %>% filter(MacroEnv == "Temperate"))
write.csv(var_table, "var_Temperate.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("variance_Temperate.png", mp, width = 6, height = 4)Mediterranean
# Save variance table
var_table <- myGxE(dd %>% filter(MacroEnv == "Mediterranean"))
write.csv(var_table, "var_Mediterranean.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("variance_Mediterranean.png", mp, width = 6, height = 4)